#ifndef INITIAL_DATA_H
#define INITIAL_DATA_H

#include "maths.h"

#include <string>

//#undef USE_RANDOM
#define USE_RANDOM
#ifdef  USE_RANDOM
//#include <boost/chrono.hpp>
#include <boost/random/mersenne_twister.hpp>
//#include <boost/random/normal_distribution.hpp>
#include <boost/random/uniform_real_distribution.hpp>

#include <chrono>
//#include <random>
#endif


/*** Base class ***/
class InitialData {
public:
  virtual real u_0(real x, real y) = 0;
  virtual real v_0(real x, real y) = 0;
  virtual real g_p_x_0(real x, real y) = 0;
  virtual real g_p_y_0(real x, real y) = 0;
  
  std::string name;
};

/*** Random generated shear layer (interface) ***/
#ifdef          USE_RANDOM
#undef          ALPHA
#define         ALPHA               0.01
#undef          NUM_ITER
#define         NUM_ITER            20
class RngInterface
: public InitialData {
public:
  RngInterface() {
    
    //std::uniform_real_distribution<double> distribution(0.0,1.0);
    //std::normal_distribution<double> distribution1(1.0,1.0);
    //std::normal_distribution<double> distribution2(0.0,1.0);
    //auto rng = std::bind(distribution, generator);
    
    
    //unsigned seed = boost::chrono::system_clock::now().time_since_epoch().count();
    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    boost::random::mt19937 generator(seed);
    //boost::random::normal_distribution<double> distribution1(1.0,1.0);
    //boost::random::normal_distribution<double> distribution2(0.0,1.0);
    
    
    boost::random::uniform_real_distribution<double> distribution1(0.0, 1.0);
    boost::random::uniform_real_distribution<double> distribution2(0.0, 1.0);
    
    
    real sum = 0;
    for(int s = 0; s < NUM_ITER; s++) {
      a_k[s] = distribution1(generator);
      sum += a_k[s];
      b_k[s] = 2 * PI * distribution2(generator);
    };
    for(int s = 0; s < NUM_ITER; s++) {
      a_k[s] /= sum;
    }
  }
  
  real u_0(real x, real y) { return (y <= atx(x) + 0.25 || y >= atx(x) + 0.75 ? 0 : 1); };
  real v_0(real x, real y) { return 0; };
  real g_p_x_0(real x, real y) { return 0.; };
  real g_p_y_0(real x, real y) { return 0.; };
  
private:
  real atx(real x) {
    real y = 0;
    for(int s = 0; s < NUM_ITER; s++) {
      y += ALPHA * a_k[s] * sin( b_k[s] + 2 * PI * s * x);
    }
    return y;
  }
  
  real a_k[NUM_ITER];
  real b_k[NUM_ITER];
};

/*** Random generated shear layer ***/
#undef          ALPHA
#define         ALPHA               0.01
#undef          NUM_ITER
#define         NUM_ITER            20
class RngShear
: public InitialData {
public:
  RngShear() {
    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    boost::random::mt19937 generator(seed);
    
    boost::random::uniform_real_distribution<double> distribution1(0.0, 1.0);
    boost::random::uniform_real_distribution<double> distribution2(0.0, 1.0);
    
    real sum_u = 0, sum_v = 0;
    for(int s = 0; s < NUM_ITER; s++) {
      a_u_k[s] = distribution1(generator);
      a_v_k[s] = distribution1(generator);
      sum_u += a_u_k[s];
      sum_v += a_v_k[s];
      b_u_k[s] = PI * distribution2(generator);
      b_v_k[s] = PI * distribution2(generator);
    }
    for(int s = 0; s < NUM_ITER; s++) {
      a_u_k[s] /= sum_u;
      a_v_k[s] /= sum_v;
    }
  }
  
  real u_0(real x, real y) {
    real prt_u = 0;
    for(int s = 0; s < NUM_ITER; s++) {
      prt_u += ALPHA * a_u_k[s] * sin( b_u_k[s] + PI * s * y);
    }
    return (y <= 0.25 || y >= 0.75 ? -1 : 1) + prt_u;
  };
  
  real v_0(real x, real y) {
    real prt_v = 0;
    for(int s = 0; s < NUM_ITER; s++) {
      prt_v += ALPHA * a_v_k[s] * sin( b_v_k[s] + PI * s * x);
    }
    return prt_v;
  };
  
  real g_p_x_0(real x, real y) { return 0.; };
  real g_p_y_0(real x, real y) { return 0.; };
  
private:
  real a_u_k[NUM_ITER];
  real b_u_k[NUM_ITER];
  real a_v_k[NUM_ITER];
  real b_v_k[NUM_ITER];
};

#endif

/*** SHEAR LAYER ***/
#undef          RHO
#define         RHO               30
#undef          DELTA
#define         DELTA             0.05
class ContShear
: public InitialData {
public:
  real u_0(real x, real y) {  return (y <= 0.5 ? tanh((y - 0.25) * RHO) : tanh((0.75 - y) * RHO)); };
  real v_0(real x, real y) { return DELTA * sin(2 * PI * x); };
  real g_p_x_0(real x, real y) { return 0.; };
  real g_p_y_0(real x, real y) { return 0.; };
};

/***  DISC SHEAR LAYER ***/
#undef          DELTA
#define         DELTA             0.05
class DiscShear
: public InitialData {
public:
  real u_0(real x, real y) { return (y <= 0.25 || y >= 0.75 ? 0 : 1); };
  real v_0(real x, real y) { return DELTA * sin(2 * PI * x); };
  real g_p_x_0(real x, real y) { return 0.; };
  real g_p_y_0(real x, real y) { return 0.; };
};

/***   STEADY ***/
class Steady
: public InitialData {
public:
  real u_0(scalar x, scalar y) { return -cos(2 * PI * x) * sin(2 * PI * y); };
  real v_0(scalar x, scalar y) { return sin(2 * PI * x) * cos(2 * PI * y); };
  real g_p_x_0(scalar x, scalar y) { return PI * sin(4 * PI * x); };
  real g_p_y_0(scalar x, scalar y) { return PI * sin(4 * PI * y); };
};

/***   VORTEX ***/
class Vortex
: public InitialData {
public:
  real u_0(scalar x, scalar y) { return -2 * sin(PI * x) * sin(PI * x) * sin(PI * y) * cos(PI * y) ; };
  real v_0(scalar x, scalar y) { return 2 * sin(PI * x) * cos(PI * x) * sin(PI * y) * sin(PI * y); };
  real g_p_x_0(scalar x, scalar y) { return 0.; };
  real g_p_y_0(scalar x, scalar y) { return 0.; };
};

#endif // INITIAL_DATA_H
